clear
clear matrix
clear mata
capture log close
#delimit ;

set mem 10000m;
set more off;

*log using ROIO22Sept15_ATCIneffectiveCourt, replace;

**************************************************************;
*  Replication Files to Accompany Ritter & Conrad ROIO 2015  *;
**************************************************************;

*The ATC figure created in this file shows the effect of the ICCPR when courts are ineffective;

use ROIOData22Sept15.dta;
sort cowcode year month;

**************************************************************;

*Generate job security variables;
gen jobsecurity_cheibub=1-haz_weibull_cheibub;
corr jobsecurity_cheibub haz_weibull_cheibub;
gen jobsecurity_full=1-haz_weibull_FULL;
corr jobsecurity_full haz_weibull_FULL;
gen jobsecurity_irregular=1-haz_weibull_irregularexit;
corr jobsecurity_irregular haz_weibull_irregularexit;

*Generate interaction term;
gen JEJS=linzerstatonJI*jobsecurity_cheibub;

*Generate dissent and repression variables;
gen dissentdummy=.;
replace dissentdummy=1 if dissct>=1;
replace dissentdummy=0 if dissct==0;
gen repressdummy=.;
replace repressdummy=0 if repct==0;
replace repressdummy=1 if repct>=1;

**************************************************************;

*Run treatment model;
capture program drop biselect;
program define biselect;
    args lnf theta1 theta2 theta3 zrhoS zrhoN;

	tempvar rhoN rhoS;
	qui gen double `rhoN' = (exp(2*`zrhoN')-1) / (exp(2*`zrhoN')+1); 
	qui gen double `rhoS' = (exp(2*`zrhoS')-1) / (exp(2*`zrhoS')+1);

    quietly replace `lnf' = ln(($ML_y1==1)*binorm(`theta1',$ML_y2*`theta2', $ML_y2*`rhoS')+($ML_y1==0)*binorm(-`theta1',$ML_y3*`theta3', -$ML_y3*`rhoN'));
          
end;

set matsize 150;
capture drop y1;
capture drop y2;
capture drop y3;
capture drop y1n;

gen y1=iccpr_rat;
gen y1n=-y1;
gen y2=2* dissentdummy -1;
gen y3=y2;

local ivar1 repressdummy linzerstatonJI jobsecurity_cheibub JEJS IO;
local ivar2 repressdummy linzerstatonJI jobsecurity_cheibub JEJS;
local ivar3 repressdummy linzerstatonJI jobsecurity_cheibub JEJS;

prob y1 `ivar1';
local ll1=e(ll);
qui matrix i1 = e(b);
qui matrix coleq i1 = eq1:;

prob dissentdummy `ivar2'; 
local ll2=e(ll);
qui matrix i2 = e(b);
qui matrix coleq i2 = eq2:;

prob dissentdummy `ivar3' if y1==0;
local ll3=e(ll);
qui matrix i3 = e(b);
qui matrix coleq i3 = eq3:;
 
heckprob dissentdummy `ivar2', select (y1 = `ivar1'); 
qui matrix i4 = e(b);
qui matrix rS=i4[1,"athrho:_cons"];
qui matrix coleq rS= eq4:;

heckprob dissentdummy `ivar3', select (y1n = `ivar1');
qui matrix i5 = e(b);
qui matrix rN=-i5[1,"athrho:_cons"];
qui matrix coleq rN= eq5:;

qui matrix initvec = i1,i2,i3,rS,rN;

ml model lf biselect (y1 = `ivar1') (y2 = `ivar2') (y3 = `ivar3') () ();
ml search;
ml maximize, difficult;

est store est1;

*The coefficient on rhos;
disp (exp(2*[eq4]_cons)-1) / (exp(2*[eq4]_cons)+1);

*The coefficient on rhon;
disp (exp(2*[eq5]_cons)-1) / (exp(2*[eq5]_cons)+1);

*************************************************************;

*Calculating predicted probabilities;
	preserve;
	set seed 10101;
	drawnorm CR_b1-CR_b18, n(10000) means(e(b)) cov(e(V)) clear;
	save simulated_betas, replace;
	restore;
	merge using simulated_betas;
	drop _merge;

	postutil clear;
	postfile mypost prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi 
            using sim , replace;
            noisily display "start";
            
	local a=0.9 ;
	while `a' <= 1.0 { ;

    {;
    
	scalar h_repressdummy=0.3342136;
	scalar h_linzerstatonJIlow=0.0191;
	scalar h_IO=3.232824;
	scalar h_Constant=1;
	
gen sel = 
	CR_b1* h_repressdummy + CR_b2* h_linzerstatonJIlow + CR_b3*`a' + CR_b4* h_linzerstatonJIlow*`a' + CR_b5*h_IO + CR_b6*h_Constant;
sum sel;

gen psel = normprob(sel);
sum psel;
centile psel, centile(2.5,97.5);

gen outcome1a =
	CR_b12* h_repressdummy + CR_b13* h_linzerstatonJIlow + CR_b14*`a' + CR_b15* h_linzerstatonJIlow*`a' + 
	CR_b16*h_Constant;
gen prob0 = normprob(outcome1a)*(1-psel);
sum prob0;
centile prob0, centile(2.5,97.5);

gen outcome1b =
	CR_b7* h_repressdummy + CR_b8* h_linzerstatonJIlow + CR_b9*`a' + CR_b10* h_linzerstatonJIlow*`a' + 
	CR_b11*h_Constant;
gen prob1 = normprob(outcome1b)*(1-psel);
sum prob1;
centile prob1, centile(2.5,97.5);

gen diff=prob1-prob0;
    
    egen probhat0=mean(prob0);
    egen probhat1=mean(prob1);
    egen diffhat=mean(diff);
    
    tempname prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi ;  

    _pctile prob0, p(2.5,97.5) ;
    scalar `lo0' = r(r1);
    scalar `hi0' = r(r2);  
    
    _pctile prob1, p(2.5,97.5);
    scalar `lo1'= r(r1);
    scalar `hi1'= r(r2);  
    
    _pctile diff, p(2.5,97.5);
    scalar `diff_lo'= r(r1);
    scalar `diff_hi'= r(r2);  

    scalar `prob_hat0'=probhat0;
    scalar `prob_hat1'=probhat1;
    scalar `diff_hat'=diffhat;
    
    post mypost (`prob_hat0') (`lo0') (`hi0') (`prob_hat1') (`lo1') (`hi1') 
                (`diff_hat') (`diff_lo') (`diff_hi') ;
    };      
    
    drop sel outcome1a outcome1b psel prob0 prob1 diff probhat0 probhat1 diffhat ;
    
    local a=`a'+ 0.01 ;
    
    display "." _c;
    
    } ;

	postclose mypost;                                  

	use sim, clear;

	gen MV = (_n-1)/100;

	graph twoway line diff_hat MV, clwidth(medium) clcolor(blue) clcolor(black) clpattern(solid)
        ||   line diff_lo  MV, clpattern(dash) clwidth(thin) clcolor(black)
        ||   line diff_hi  MV, clpattern(dash) clwidth(thin) clcolor(black)
        ||   ,   
            xlabel(, labsize(4)) 
            ylabel(-0.1 -0.05 0 0.05 0.1 0.15 0.2, labsize(4))
            yscale(noline)
            xscale(noline)
            yline(0)
            legend(off)
            title("", size(3))
            subtitle("" "" " ", size(3))
            xtitle("Executive Job Security", size(3)   )
            ytitle("Change in Pr(Mobilization)", size(3))
            xsca(titlegap(2))
            ysca(titlegap(2))
            scheme(s2mono) graphregion(fcolor(white));
			gr_edit xaxis1.reset_rule 3, tickset(major) ruletype(suggest);
			gr_edit .xaxis1.major.num_rule_ticks = 3;
			gr_edit xaxis1.edit_tick 2 0.05 `"0.95"', tickset(major);
			gr_edit xaxis1.major.num_rule_ticks = 2;
			gr_edit xaxis1.edit_tick 3 0.05 `"0.95"', tickset(major);
			gr_edit xaxis1.major.num_rule_ticks = 2;
			gr_edit xaxis1.edit_tick 2 0.1 `"1.0"', tickset(major);
			gr_edit xaxis1.major.num_rule_ticks = 1;
			gr_edit xaxis1.edit_tick 1 0 `"0.9"', tickset(major);
			gr_edit xaxis1.reset_rule 3, tickset(major) ruletype(suggest);

graph export ICCPRAT_IneffectiveCourt.eps, replace;

*************************************************************;

exit;
